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Abstract 

We present the first exact analytic result for all partonic channels contributing to the total cross 
section for the production of a pair of heavy flavors in hadronic collisions in NLO QCD. Our 
calculation is an essential step in the derivation of the top quark pair production cross section 
at NNLO in QCD, which is a cornerstone of the precision LHC program. Our results uncover 
the analytical structures behind observables with heavy flavors at higher orders. They also reveal 
surprising and non-trivial implications for kinematics close to partonic threshold. 



1 Introduction 



Thanks to the Tevatron, the production of heavy flavors at hadron coUiders has nowadays become 
a very dynamic field of research. The imminent physics expected from the Large Hadron Col- 
lider (LHC) places even stronger emphasis on this class of processes. The very large production 
rates for both top and bottom at the LHC will alow studies of heavy flavors with precision so far 
seriously considered only for processes with light flavors. Accurate knowledge of top produc- 
tion is fundamental in Higgs and most Beyond the Standard Model as well as Standard Model 
studies [1-5]. 

These very exciting experimental developments have to be confronted by theory. At present, 
the next-to-leading order (NLO) corrections in the strong coupling have been accounted for. That 
includes the fully inclusive production rate, the one-particle inclusive cross section, spin correla- 
tions as well as the fully exclusive production and charge asymmetry [6-15]. At the NLO accuracy 
level several associated production processes have also been studied [16-24]. 

Various sources of uncertainty affect the theoretical predictions for heavy flavor hadroproduc- 
tion. The most important one is due to unknown higher order corrections; these are NNLO effects 
in the strong coupling or electroweak effects. Their reduction mandates higher order calculations. 

Partial NNLO corrections to the total cross section have been presented in the literature [25- 
27]. They originate from a fixed order truncation of the all-order exponentiation of soft gluon 
logarithms [28-32]. We would like to emphasize, however, that such an approach to the fixed 
order cross section does not represent a controllable approximation in the sense of being derived 
from a consistent expansion in a small parameter. A priori, such partial corrections cannot be 
expected to approximate well the exact fixed-order result; an explicit example can be found in 
Fig. 3 in Ref. [33] as well as in Section |4]below. A careful investigation of this point is beyond the 
scope of the present article. It will be provided elsewhere [34]. 

In a recent analysis, the authors of Ref. [35] point out that bound-state effects may have been 
underestimated in past studies of kinematics very close to threshold. They argue that this effect 
can have a significant impact on both near-threshold differential distributions and top mass mea- 
surements. 

Another source of uncertainty in top production are initial state non-perturbative effects con- 
tained in the parton distribution functions. These effects have been discussed extensively in the 
literature; for a comparative study of the various pdf sets in top quark pair production see for ex- 
ample Ref. [26,32]. The importance of top quark pair production in constraining pdf sets and 
other LHC observables has been studied in Ref. [36]. Overall, pdf -related effects in the considered 
process are expected to be in the 3-4 percent range and are smaller than the missing higher order 
effects of interest to us in this work. 

The present paper represents a needed step in the program [37-43] for the systematic derivation 
of the NNLO QCD corrections to the total heavy flavor production cross section. The phenomeno- 
logical importance of the total pair-production cross section was established in the original pa- 
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pers [8], where it was demonstrated that the effect of the NLO corrections is mainly on the overall 
normalizations and not on the shapes of differential distributions (however see Ref. [35] as men- 
tioned above). It is natural to expect that this feature will persist at even higher orders, especially 
for distributions not subjected to too strong cuts. The reason for this expectation can be traced 
back to the non-trivial kinematics of the produced heavy quark pair starting from the leading order 
approximation. 

Specifically, in this paper we derive the first exact analytic expression for the NLO QCD cor- 
rections to the total inclusive cross section for heavy quark pair production at hadron colliders. The 
results obtained are thus new, since previously only a numerically extracted fit was available. We 
comment on the quality of this numerical approximation in Sec. IH Secondly, we aim at tuning a 
calculational approach that is capable of tackling the derivation of the NNLO corrections for this 
observable. 

The calculation of the top quark pair production cross section at higher orders presents several 
significant challenges. First, the analytic structure of the total cross section is very complicated - 
and different from the one encountered so far in multiloop massless calculations. Indeed, radicals 
start to show up even at leading order. Our approach in addressing this problem is discussed in 
Section[3l Our findings (see sections. [4I5I6I) not only confirm the expectations, but also demonstrate 
that a priori unforeseen analytic structures arise starting from NLO. Second, at NNLO one also 
has to tackle the problem of the large number of three-loop cut diagrams. Third, in the case of 
top quark pair production the value of the top mass is not small compared to the typical partonic 
center of mass energy. That in turn means that one has to obtain results, which are exact in the 
mass. In consequence, approaches based solely on an expansion in powers of the mass of the heavy 
flavor are not likely to be sufficient to cover the whole kinematic range. We illustrate this point by 
analyzing the newly derived exact NLO result in Sec. IH 

Finally, we would like to address the question of why we prefer an analytic approach and do 
not pursue a purely numerical one from the very beginning. One reason is that the NLO result, 
including terms subleading in the regularization parameter, is a prerequisite for the NNLO calcu- 
lation. It is clearly desirable to have complete control over the lower order result. In fact, as we 
show in our following discussion, effects due to uncertainties inherent in numerical calculations 
can be sizable enough to require additional consideration. Our second motivation for pursuing 
analytical methods is based on the observation that the cross section we compute essentially rep- 
resents a one-variable problem. Therefore, it seems reasonable to expect that at least most of the 
calculation should be feasible in terms of known functions and by using automatic methods. Our 
expectations have indeed turned out to be true. As a result we have gained deep insight into the 
analytic structure of cross sections as complex valued functions in the massive case. 

The paper is organized as follows: in the next section we introduce some basic notations and 
definitions. In Sec. [3] we introduce our calculational method. Section |4] contains our main result, 
the analytic expressions for all partonic channels contributing at NLO. Sections [5] and [6l present 
detailed discussions of the analytic features of our results. Finally, in order to be as self-contained 
as possible, we collected various known formulae in the Appendix. 
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2 Setup 



We calculate the NLO corrections to the total cross section for the process 

p+p{p) ^ e+(2+^- (1) 

At partonic level, this involves the following reactions 

q+q 

g+g ^ 2+2+^, 

g + q ^ Q + Q + X. (2) 
Following the established notations, we denote the cross section for the reaction Eq. ([1]) as 

o(5,m2) = / dxidx20ij{s,n?,i/)fi{xi,i/)fj{x2,i/) , (3) 

where s = X1X2S is the partonic center of mass energy, S is the corresponding hadronic invariant 
and m is the pole mass of the heavy quark. The physical region is defined by 5 > 4m^ > 0. The 
scale /J denotes both the renormalization and factorization scales as appropriate. For simplicity, 
we do not distinguish between the two. The functions fi are the parton distribution functions for 
parton i in the (anti)proton. 

The partonic cross section can be written as an expansion in the strong coupling. We will keep 
the traditional notation, where the cross sections are written in terms of the following dimensionless 
functions 

2 2^ «s(^^) r 



Oijis^y) = ^^fij i—^h (4) 



where 

.2 ,,2 



f-'C^^^ (t) — (4" (t) (t)'-(^)) - ■ 

The partonic cross section a contains three implicit scheme dependencies. First, it depends on 
the definition of the strong coupling a^. We work in terms of the standard MS coupling defined 
through 



(6) 



e V27r 

where 5e = (47c)^exp(— £Ye) and |3o = — ^Tftij: is the first term of the QCD |3-function (known 
by now up to the four-loop level [44,45]). The color factors in an SU(A/^)-gauge theory are Ca ~ N, 
Cp = (A^2 _ i)/(2N) and Tp = 1 /2. Throughout this paper stands for the total number of flavors, 
which is the sum of light and rij^ = 1 heavy quarks. Switching between definitions of the coupling 
with + 1 (as adopted in this paper) and active flavors is done with the standard decoupling 
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relations (also known up to the four-loop level [46,47], see Appendix lAl for the expression needed 
in our case). 

Second, the partonic cross section d depends on the definition of the mass m. As we mentioned 
above, we work in terms of the usual pole mass [48,49]. The third implicit dependence is on 
the scheme used to factor out the coUinear (mass) singularities related to initial state coUinear 
radiation. We adopt the standard MS factorization scheme defined through 



I 

kl 



® rfc,-(x,£) 



r/,(^,e) 



(p), 



(7) 



where p = Am^ j s with < p < 1 . The variable x, p < x < 1 , stands for the corresponding (dummy) 
convolution variable. The functions Oyt/ are the "bare" (i.e. coUinearly un-renormalized) hard 
scattering cross sections and the coUinear counterterms F are given by 



r,;(x,£)=5,,5(l-x)- 




+ o(a2) 



(8) 



where (x) are the {n + l)-loop (space-like) QCD splitting functions [50-57]. The superscript 
in P^") and tts emphasizes that these quantities are evaluated in a scheme with n^-flavors. 



3 Parton level evaluation 

In this section we describe the evaluation of the partonic cross section o,y(p, £) appearing in Eq. O. 
Following the method introduced in [58] we perform simultaneously all loop and real radiation 
integrations. With the help of the IBP identities [59,60] which we solve with the Laporta algorithm 
[61], the whole problem is mapped into a total of 37 master integrals. 

The needed partonic cross sections are non-trivial functions of a single dimensionless variable 
As was established in the literature, the analytic structure of the result becomes 
rather complicated. In particular, even at LO radicals like P = ^1 — p appear which suggests that 
the application of established analytic methods in this problem may be too cumbersome. To rectify 
the situation we perform the following change of variables 



s (1+x 



i2 



1_ /l_4mi 

X = (9) 

i+Vi-¥ 

This change of variables allows one to express the results for the partonic cross sections almost 
entirely in terms of harmonic polylogarithms (HPL) [62] of the variable x as well as rational func- 
tions of X and IrLx. Such analytic structures have been widely explored and can be dealt with in a 



'The dependence on the renormalization/factorization scale /j^/m^ requires only calculations of one order lower 
and are thus considered 'trivial'. However, see the discussion at the end of this section. 
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straightforward manner. A particularly welcoming feature of this approach is that one can achieve 
a large degree of automatization in the solving of the differential equations for the masters with the 
help of existing software [63,64]. 

The reason why this change of variables is helpful can be seen as follows. The physical cross 
section for the production of two heavy particles in the final state has two physical singularities: 
the points = representing the high energy (or massless) limit, and the kinematic threshold 
for heavy pair production s = Anp-. The third singularity, which is not in the physical region, is 
\m'^ / s\ oo. The transformation ^ maps these three singularities into the points x = (0, 1, — 1) 
which are precisely the singular points of the differential equations defining the set of HPL's. 

As we will see in the following, however, the differential equations for a few of the master 
integrals posses singularities other than the three described above. These are the points 5 = m^, 5 = 
—nP'.s = —Am} and s = —ISir}. They are all outside the physical region < 4m^/s < 1. To the 
best of our knowledge, this is the first calculation of an observable with massive flavors where 
such pseudo-threshold-like singular points have been encountered. At NLO they appear only in 
the g + g ^ QQ + X subprocess starting from order □(e'^). We discuss them in detail in SecjS] 

To completely fix the solutions for the master integrals one has to supply suitable boundary 
conditions. The behavior of the solution around any one of the singularities x = (0, ±1) can be 
used. In the cases where integrals form a system of two or three differential equations it turns out 
that one needs to supply information about the behavior of the solution at more than one singular 
point. It is also worth mentioning that the approach used in Ref. [65, 66] to fix the boundary 
conditions by simply requiring regularity of the solution cannot be applied to the problem at hand 
since the masters are typically singular at all singular points. A detailed discussion of the most 
complicated master we encountered in this calculation is given in Sec. |6l Here, let us only note 
that while some of the boundaries required a direct calculation by Mellin-Barnes methods [67,68], 
where we used the MB package [69], many could be simply obtained by requiring the vanishing of 
the integral at threshold. 

The evaluation of the convolutions in Eq. (|7]) also deserves a comment. A complication arises 
from the fact that while the dummy variable x is the 'natural' variable for the splitting functions 
entering the coUinear counterterms F, the same variable is not the natural one for the partonic cross 
sections a. To write the latter in their natural variables one has to switch from the variable p to the 
conformal variable x defined in Eq. 

4 Results 

At this point, we are ready to present our main result, namely the analytic expressions for the 
partonic cross sections at next-to-leading order. Our calculation has been performed with full 
account of the dependence on the number of colors A^. Because of the substantial length of the 
formulae, however, we will present here only the result for A'^ = 3. The complete results, valid for 
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any SU (N) color gauge group are available in electronic form in [70]. 
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Figure 1: Error implied by the numerical integration of the cross section [6] in the quark-anti- 
quark production channel against the exact result. The functions are defined in Eqns. U0\25^ . 
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Figure 2: Error implied by the numerical integration of the cross section [6] in the quark-gluon 
production channel against the exact result. The functions are defined in Eqns. Ul\25^ . 



+ 
+ 



8 



+ ■ 



32 



16 



9{x+l) 2,{x+\Y 9{x+\f 9(jc+l)4 

16 16 64 32 

+ 



9{x+l) 3(x+l)2 9{x+lf 9{x+lY 



H(0,0,0,.j;) 
H(0,0,l,x) , 



(11) 



+ 



+ 



+ 



+ 



13121 



+ ■ 



38761 



1506193 



+ ■ 



312977 



237377 



21 



V1920(.x+1) 1152(a;+ 1)2 8640(.x+1)3 \AAO{x+\f 3600(a;+1)5 {x+\f 
5 7291 9 



+ 



{x+iy 5400 64(;c-l) 



x + 2 



+ 



11227 



81 



13999 



27 



202939 



n{x^+x+\) 32(a:- 1) 576(a;+1) 32(a;-1)2 6A{x+\Y \6{x-\f 288(.x+l)3 



80471 



+ 



80471 



96{x+lY 240(;t;+l)5 180 y 



181\ , 
H(-l,.x 

1 OA I ^ ^ > 



x-l 



24665 



+ ■ 



5x 



■ + ■ 



301319 



881899 1009037 

■ + 



18(x2+x+l) 576(x+l) 18(x2 + 3x+l) 576(x+l)2 576(x+l)3 576(x+l)4 

24 6 181 9 9 



1173109 2311 
+ ■ 



+ ■ 



181 



+ ■ 



+ 



1440(x+l)5 18(x+l)6 {x+iy {x+lf 60 8(x-l) 4(x-l)2 
27 \ 

H(0,x) 



32(x- 1)3 
181 260 



+ ■ 



3784 



12898 15563 

■ + ■ 



31126 



90 9(x+l) 9(x+l)2 9(x+l)3 9(x+l)4 45(x+l)5 
5(3x+l) 27 10165 5(llx + 2) 



H(l,x) 



+ ■ 



+ ■ 



135 



+ 



108(x2 + 3x+l)^ 128(x-l) 3456(x+l) 216(x2 + 3x+l) 128(x-l)2 
26279 81 17479 81 34091 5999 



+ ■ 



3456(x+l)2 32(x-l)3 1728(x+l)3 64(x-l)4 864(x+l)4 192(x+l)5 



8 



10049 ^ ^ 2 



1728(a:+1)6 4{x+l) 
( x+\ 63 llx-l 5(llx + 2) 315 5(3;c+l) 



^6(;c2+;c+lf 32(^-1) 12(^2+^+i) i8(;c2 + 3^+l) 32(^-1)2 9(;,2 + 3^+i) 
189 189 607 14245 5263 6533 221 



2 



8(;t:-l)3 16(;c-l)4 288(;c+l) 288(x+l)2 36(x+l)3 72(;c+l)4 6(x+l)5 

593 \ , 

( 20 220 310 124 \ 

( x+l 9 llx-2 45 27 27 833 

+ — ^-TTT TT + TTTT^^ ^ + 7TT + :TT + ' 



^6(;c2+x+lf 16(x-l) 12(x2+x+l) 16(x-l)2 4(x-l)3 8(x-l)4 48(x+l) 

2429 11425 16555 2725 _ 1733 \ 

^72(x+ 1)2 + 144(x+ 1)3 " 72{x+lY ^ 16{x+ 1)5 " 48(x+ 1)6 J ^ '""^ 

2-llx 45 10769 5(llx + 2) 225 23599 



12(^2 +X+1) 64(x-l) 576(x+l) 36(^2 + 3x+l) 64(.x-l)2 576(x+l)2 
-x-l 5(3x+l) 135 19513 135 32605 



6(;c2+x+l)' 18(x2 + 3x+l)' 16(x-l)3 288(x+l)3 32(x-l)4 144(x+l)4 



16219 2693 

+ 



9 \ 



96{x+iy 96{x + 
. 361 1193 2354 1291 1331 773 \ , 

^ \9ix+ 1) " 18(x+ 1)2 " 9(x+ 1)3 + 2(x+ 1)4 " 3(x+ 1)5 + 9(x+ 1)6 ) ^ ' '""^ 
14 20 220 310 124 \ 

751 1069 983 4142 545 13 ^ 

V36(x+ 1) ^ 36(x+ 1)2 " 3(x+ 1)3 + 9(x+ 1)4 " 3{x+ 1)5 " 9(x+ 1)6 J ^ ' '""^ 
28 40 440 620 248 \ 

3433 5935 25709 7405 385 2897 20 



576(x+l) 144(x+l)2 288(x+l)3 144(x+l)4 96(x+l)5 288(x+l)6 9(x+l) 
571 81 567 405 81 \ 

^192(x-l) ^ 32(x-l)2 + 16(x-l)3 ^ 8(x-l)4 + 4(x-l)5 j 

/ 997 2839 4771 5221 5 5 \ ^ . . 

^ V432(x+ 1) ~ 216(x+ 1)2 + 216(x+ 1)3 ~ 432(x+ 1)4 + 4(x+ 1)5 " 12(x+ 1)6 ) ^ "("^'"^^ 
8375 61711 56633 4195 1625 269 

+ 1- ^^.... ..^ + .^.^. ..^. +■ 



7 



27648(x+l) 6912(x+l)2 3456(x+l)3 576(x+l)4 864(x+l)5 216(x+l)6 
5 1693 81 567 405 81 \ 2TT/f^ 

^27(x+l)^ ^9216(x-l) ^512(x-l)2 + 256(x-l)3 + 128(x- 1)4 + 64(x- l)5y) ^ ^ '""^ 

/ 43 449 212 164 36 ^ 1 1 ^ 

^ V 6(x+ 1) ^ 18(x+ 1)2 " 9(x+ 1)3 " 9(x+ 1)4 + (x+ 1)5 " (x+ 1)6 J ^" ' '""^ 

/ 545 154 1429 2059 21 ^ 1 1 ^ 

^ V36(x+ 1) " 3(x+ 1)2 + 18(x+ 1)3 " 36(x+ 1)4 + (x+ 1)5 " (x+ 1)6 J 
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/ 1187 2581 1439 2009 57 19 \ tt/ i n n 

^ V 72(.T+ 1) + 36(;c+ 1)2 " 12(;c+ 1)3 + 2A{x+ If ~ 2{x+ 1)5 ^ 2(;c+ 1)^ J ^^"^'"'"'^^ 
/ 25 117 184 92 \ , 

+ 1-7 + 7 H(-1,0, l.x) 

V ^+1 (^+1)2 [x+iy^ {X+\Y J ^ ' ' ' ^ 

/ 23447 3713 11467 7213 1345 1303 40 



V2304(;c+ 1) 576(a;+1)2 288(;c+1)3 144(a:+1)4 36(a;+1)5 72(;c+1)6 9{x+\) 

1315 189 1323 945 189 \ 

^768(x- 1) + 128(x- 1)2 + 64(x- 1)3 + 32(;c- 1)4 + 16(;c- 1)5 ) H(0,-1,0,-^) 

/ 8 24 56 8 24 ^ ^ 1 1 ^ 

4561 10313 9313 1081 167 3 55 
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1152(x+l) 288(x+l)2 144(x+l)3 36(x+l)4 72(x+l)5 8(x+l)6 128(x-l) 

27 189 135 27 \ ^ 

H(0,0,-l,x) 



64(x-l)2 32(x-l)3 16(x-l)4 8(x-l)5 
/ 14875 18817 54893 5851 241 427 10 

+ L.^.. .,.. + .^^, ..^. + 



V4608(x+1) 384(x+l)2 576(x+l)3 96(x+l)4 16(x+l)5 72(x+l)6 9(x+l)7 
1001 135 945 675 135 



1536(x-l) 256(x-l)2 128(x-l)3 64(x-l)4 32(x-l)5 



H(0,0,0,x) 



43 377 3085 392 172 248 40 



72(x+l) 4(x+l)2 18(x+l)3 9(x+l)4 3(x+l)5 9(x+l)6 9(x+l)7 
^ H(0,0,l,x) 



24(x- 1) 



24 56 8 24 ^ ^ m 

1735 289 3107 173 76 38 20 

+ \- ..., .... + ^.. .... + ^, .... +■ 



144(x+l) 8(x+l)2 36(x+l)3 9(x+l)4 3(x+l)5 9(x+l)6 9(x+l)7 

/ 1 5 5_ 1 \ 

^ V 12(x+ 1) ^ 12(x+ 1)3 " 6(x+ 1)4 + 4(x+ 1)5 " 4(x+ 1)6 ) ^^""^ 

( 5 5 5 10 5 ^ \ ( \ 

^ V 108(x+ 1) ^ 54(x+ 1)2 + 108(x+ 1)3 " 27(x+ 1)4 + 12(x+ 1)5 " 36(x+ 1)^ ) ^^^""^ 
4 4 16 112 16 16 128 



+ : 



27 27(x+l) 9(x+l)2 27(x+l)3 27(x+l)4 3(x+l)5 27(x+l)6 
32 \ F3(x) 



27(x+l)V Vx^ + 6x+l 
+ F4(x) 

r x^{-2x+{x+l)U{0,x)+2) \ 



Most of the results above are expressed in terms of the standard harmonic polylogarithms 
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1.5 




Figure 3: Error implied by the numerical integration of the cross section [6] in the gluon-gluon 
production channel against the exact result. The functions are defined in Eqns. di2l25D . The plot 
on the right contains also six terms of the threshold (long-dashed line) and the high energy (short- 
dashed line) expansions. 

(HPL) [62] of weight up to three and of the simple argument x, which implies that they can also 
be rewritten in terms of the usual polylogarithms of more complicated arguments but same weight. 
As we already mentioned in Section[3]and further elaborate in Section[5l the presence of additional 
singularities beyond those naively expected makes it impossible to integrate the gluon-gluon cross 
section in terms of HPL's alone. Indeed, four additional functions are need. The first three of them 
are chosen as follows 

[' (2z+l)(H(-l,0,z)+H(0,-l,z)-H(0,0,z)) 
Jx 2 z^+z+1 



1 (2z + 3) ( 12 H (- 1 , 0, z) - 6 H (0, 0, z) + 7t^) 
4(z2 + 3z+r 



= -/ Az'-^^ — \,2:i:.. ^-----^ (14) 



^f' 5(z-l)(l2H(-l,0,z)-6H(0,0,z)+7t2) 
¥^(x) = + dz ^ , , -. (15) 

It is clear from these expressions that both Fi and F2 can also be expressed through the usual poly- 
logarithms of up to weight three. Since the resulting expressions are very lengthy and cumbersome, 
yet obtainable by automated software, we prefer to stick to the integral representations which are 
easier to manipulate. Unlike Fi 2, however, the presence of the square root in the integrand of F3 
makes it impossible to give a closed form result in terms of standard special functions in that case. 

Besides being able to provide numerical values for the partonic cross sections at any given 
point, we would also like to perform threshold and high energy expansions. As far as HPL's 
are concerned this can easily be done with the help of [63]. As it turns out, Fi, F2 and F3 can 
be expanded with little additional effort by simply integrating the expansions of the integrands 
term-by-term. The missing boundaries are trivial at threshold, since all integrals vanish there by 
their very definition. In the high energy limit, the asymptotic behavior necessary to complete this 
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procedure is 



/m2\ 13^3 (n?\ 

Fi — = ^ + — , (16) 



s ) 24 V 

/m2\ 21^3 

F2 — = ^ + — , (17) 



s ) 10 V * 

F3 — = -ologM— +^log — +-ri^ + o — . (18) 



The origin of the last function, F4, is discussed at length in Section [6l We choose to write its 
integral representation in terms of the original mass variable p = Am^ / s, not x, (recall that the two 
are related through Eq. ^) because this simplifies the arguments of the special functions. We have 

F4(x)= /''dx/4(p, X), (19) 



where 



45p 



^4(p,x) = t;^;— log 



l-v^\ /((p^ + 1)K(v/=?p)-(p-1)E(v^))k(-^ 



3271X ''Vi+v^r^y I V4xTT 

( (-4p2 + 3p + 1) E (^) + (3p2 - 3p - 2) K (^) ) K (^) 

v/4pTT 



(20) 



The integrand is expressed through the complete elliptic function of the first, K, and second, E, 
kind defined as 

- r^vi-zvi-.v '''' 



m = /'dz^^^iJ^, (22) 
Jo VI -z^ 

Obtaining the asymptotic behavior of F4 in the two interesting limits is a substantially more 
complicated task than in the case of Fi . . .F3. This is best done with the help of differential equa- 
tions. Here we only reproduce the leading term of the expansions 

69 



F4(|3) = -^13^ + 0(136) , (23) 



and 



^o('^). (24, 
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Having given all of the relevant expressions, we would now like to address the question of the 
actual precision of the fitting formulae derived by numerical integration in [6]. To that end, we plot 
the error defined as 



Af p (P, = 0) = ''^ ' '-^ ^ , (25) 

where f^p'^^^ has been taken from [6], as function of the velocity, |3, over the full range from the 
threshold (3 = to the high energy limit P = 1. As can be seen in Figs. (I1I2I3I) the error exceeds 
the promised 1% only in the vicinity of a zero of the cross section correction. This is of course 
expected, since there is no particular need or phenomenological interest to find the exact position 
of the zeros only to have an improved fitting formula, and the authors of [6] have not used this 
information to obtain their approximations. Overall, as can be concluded from these plots, the 
quality of the fits in the case of fj^^^ and fgl^ is excellent. However, the error almost approaches 

the 1% bound in the case of/^]\ While some loss of precision is expected given how complicated 
a function the gluon-gluon cross section is, we will see in the following that this slightly larger 
uncertainty is of relevance to the threshold behavior of the cross section. 

The analytic results presented above make it possible to obtain exact expansions in both the 
massless and the threshold limits. As far as practical applications at colliders are concerned, the 
threshold behavior is more important. We give it here with full account of the color factors and Tf 
set to its customary value of 1/2 

4l\^) = 4^/i°^(P)((cF-^Q)|^ + 2CFlog2(8|32)-(8CF + Q)log(8|32) 
+ CF(^8-y + 31og2-21og22^+CA(^y-y-21og2 



.9 3 J 9 



- 4^/r(P) C.-A^^^U^ + 2Qlog^(8P^)- - - log(8P^) 



/21iV2-50 (l7iV2_4o)^2 (A^2_4)i„g2 
+ ^^[-1^ 24(A^2_2) 21og2 

+ CF(^-5 + y^+0(p)^ . (27) 

The expressions for the Bom level functions fj^^^ and fgg^ can be found in Appendix lAl The singu- 
lar term proportional to 1 /(3 (Coulomb enhancement) as well as the coefficients of the logarithms 
log(8|3^) (soft-gluon enhancement) have been known already from [6]. 

The (3-independent constant terms are important for the quantitative understanding of heavy- 
pair production at threshold. Numerical approximations were first derived in [6] . Recently the au- 
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thors of Ref. [35] independently extracted these terms from earlier calculations [71,72] on quarko- 
nium production at hadron colliders. The comparison of the results extracted in [35] with the one 
derived from the calculation of [6] showed that the two results were incompatible with each other, 
both for qq and gg scattering. The numerical size of the discrepancy was found to be significantly 
larger than the numerical uncertainty inherent in the results of Ref. [6]. 

With the help of our calculation this discrepancy is now resolved and the conclusions are quite 
intriguing. In the qq channel the discrepancy turned out to be due to a missing contribution in the 
quarkonium calculation, see Ref. [35] for details. The situation is quite different in the case of 
the gluon-gluon channel. In this reaction our calculation agrees with the analytic result extracted 
in Ref. [35]. Thus, the source of the discrepancy is solely due to an accidental magnification of 
the numerical uncertainty of the calculation in Ref. [6] . Indeed, the numerical counterpart of the 
^-independent constant term in Eq. (|27|) 

^ ( 21N^-5Q {\lN^-AQ)ie , (A^^-4)log2 \ / 

""^-W^- 24(iV^-2) +^V^ 21og2j+C.^-5 + - 



21 168 
is according to [6] 



nil 28371^ 151og2 ^, 2. _ 

+ — -^-61og2 2~ 34.88, (28) 



76871 



^ flf -37.25. (29) 

The difference between the two results is about 7%, even though the fitting formula for the cross 
section from Ref. [6] is more accurate than 1%, see Fig. [3] This is not entirely surprising as the 
following arguments show. Close to threshold the term proportional to this coefficient in fg^g is 
damped by a factor of |3, whereas the dominant contribution comes from logarithms of the velocity 
and Coulomb enhanced constant terms. Thus its error is less relevant for the final precision. For 
larger values of [3 other terms dominate and their errors conspire to make the result sufficiently 
precise. An analysis showing that this mechanism indeed works, i.e. changing the value of Oq by 
up to 10% leaves the total cross section within its 1% bounds, has been performed in [35]. 

This difference represents a powerful example for the need of analytic results in certain cases. 
It will be very interesting to further study how this effect propagates in the soft-gluon resummed 
cross section. Work is in progress [34]. 

Let us finally turn our attention to the high energy limit of the cross sections. It is well known 
that the ^-channel gluon exchange leads to a constant behavior of the quark-gluon and gluon-gluon 
cross sections in this regime. We derive 

The constant in Eq. (|30l) is in perfect agreement with the numerical value obtained in [6] . The 
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relation between the two cross sections as given in the second equation above has been explained 
there as well. 

The quality of deeper expansions around P = and (3 = 1 is shown in Fig.[3l With six terms in 
both cases, most of the range is covered satisfactorily. In fact, due to large numerical cancellations 
in the vicinity of P = 0, it is better to use the expansion rather than the exact formula. Let us stress 
that these cancellations are a problem only in the case of numerical evaluation. Higher precision 
in the numerics always allows to retrieve the correct result. 

5 Singularities of the differential equations 

As already mentioned in SeclH besides s = 4m^, = and s = 0, we have observed four addi- 
tional singular points of the differential equations of the master integrals. These additional singu- 
larities characterize the integrals depicted in Fig. HI which contribute exclusively to the gluon-gluon 
cross section. 

The presence of additional singularities in the differential equations may have two different im- 
plications. One possibility is that, in principle, a different choice of master integrals may remove 
these singularities if they do not have some deeper (physical) meaning. However, if these singu- 
larities persist in the amplitude, it immediately follows that it would not be possible to express the 
final result through harmonic polylogarithms of the x variable alone. 

We study these alternatives in the following. The main idea behind our analysis is that, although 
we work with cut graphs, the homogenous parts of the corresponding differential equations are the 
same as for ordinary Feynman integrals, which means that if the corresponding Feynman integral 
has a singularity then the cut graph may have it as well. 

It turns out that the singularity at s = —I6m^ occurring in the graph of Fig. |4] a) is the one 
easiest to explain. Indeed, a quick look at the Landau equations shows that this is the leading 
singularity of the Feynman integral, when all lines go on shell. As such, it can be classified as 
an anomalous pseudo-threshold, because it occurs in the Euclidean region. How to integrate the 
differential equations in this case is explained in Section [6l 

The second singular point, which is also rather easy to explain, is 5 = — 4m^. Let us take a 
closer look at Fig.|4]b). The graph is crossed in the forward direction because of the exchanged 
outgoing momenta. Its existence is clearly due to the possibility to switch the identical external 
gluons. If we treat the graph as a normal box graph with arbitrary external momenta, then it would 
have a physical singularity at t = 4m^. The forward limit of the crossed graph corresponds here 
to the Mandelstam u = 0, and not t = 0. This in turn implies that t = —s and the usual branch 
point modeled by a logarithm is mapped to —s = 4m^. Interestingly, this argument further implies 
that this singularity can occur in any crossed box graph (in the sense explained above), because 
the heavy quark line always has to form a closed loop, and therefore the singular point t — An? is 
always present. 
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a) 



s = —IQw? 




Pi 



P2 



s = — 4m^ 



P2 



Pi 



Pi 



P2 



s = m 

P2 
Pi 



Pi 



P2 



-m 



P2 



Pi 



Figure 4: Cut forward Feynman integrals whose differential equations introduce a priori unex- 
pected singularities at the listed points. Thick lines are massive, whereas thin are massless. Dashed 
lines represent cuts (cut lines are on-shell). Notice the crossed momentum flow in the b), c) and d) 
cases. 



Interestingly, the seemingly simplest graphs in our problem have the singularity that is hardest 
to explain. The singularity of Fig. |4]c) at s = m^ does not correspond to any leading singularity 
in the sense of Landau equations. If we kept the external momenta as free, then there would be 
three master integrals for this topology and mass distribution. Their differential equations do not 
have more singularities than those of the one-loop graph without the massless internal line. It is 
only when we reach the forward direction that the three integrals collapse to one and generate the 
aforementioned singularity in the process. Fortunately, Fig.|4]c) and d) are linked. In fact going 
from one to the other while forgetting the cut, can be done by the change s ^ t. Since in the 
forward direction t = —s, the singularity at s = m^ of c) implies a singularity of d) at 5 = —m^. 

A closer examination of the final result for the gg cross section Eq. (fT2l) shows that all four 
singular points described above are present in the final result in some form or another. Indeed, 
Figs. |4]a), b), c) and d) correspond to F4, F3, Fi and F2 respectively. Whereas s = —I6m~ for F4 
and s = —m^ for F2 are usual branch points of logarithmic type, s = — 4m^ is only a square root 
type branch point of F3, and s = m^ is not a branch point of Fi at all. This is exactly the opposite 
of what would happen if we had to do with an amplitude, because there, all points in the Euclidean 
region would have to be regular. 

As a final remark we would like to stress that while at next-to-leading order these additional 
singularities are only present in the gluon-gluon channel, we anticipate that they will be present 
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Figure 5: The topology and mass distribution of the most complicated Feynman integrals, together 
with the momentum distribution necessary to define the numerator of the integrand considered in 
the text. Thin lines are massless, whereas thick have the same non-vanishing mass. The dashed 
line represents a cut (cut lines are on-shell). 

also in the quark-anti-quark channel at the next-to-next-to leading order. 



6 The most difficult master integrals 

Of all 37 master integrals occurring in our problem, two stand out as particularly difficult. They 
both correspond to the topology and mass distribution depicted in Fig. [5] and contribute to the 
gluon-gluon cross section exclusively. As usual in such cases, integration-by-parts (IBP) identities 
make it possible to choose any powers of the propagators and irreducible numerators, as long as 
one of the integrals cannot be reduced to the other (independence with respect to IBP's). We make 
the following choice 

Ml = s'-'r^{l-e) [^^-^8^{k^-m^)b^{r^-m^)8+{l^) ^ (32) 



M2 = s'-'Y\l - E) / [k^ - m2) 5+ (r2 - m^) 5+ {f) , (33) 



where are Feynman propagator denominators of the lines which are not cut, pi and r are defined 
in Fig[5l whereas k and / are the momenta of the remaining cut lines also shown in Fig. [5l The 
sign of the 5- functions defines the energy flow of the cut lines, as usual. 

One can show that both M\ and M2 are finite (with or without the cut). Moreover, they enter the 
cross section with coefficients regular in the dimensional regularization parameter £, which means 
that we only need to know their value at £ = 0. In terms of the dimensionless variable p = Am^ /s, 
they satisfy the following system of differential equations 

d_m ^ _^1±1m, + ^M, + H,^ (34) 
^ dp 4p + l ^ 4p + l ^ ^ ' ^ ' 
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dM2 3p + 1 2p + 1 

p— ^ = - o.. ,, Mi + -^^M2+//2. (35) 
Jp 8(4p + l) 4p + l 

Although their exact form is irrelevant for our forthcoming discussion we reproduce here the non- 
homogenous terms H\ and H2 

5 (log (fe^) +2^1^) 

- 4^T1 ' ^^^^ 

= 4^1. (37) 
Id 

Let us rewrite the system as one second order equation for M\ . We obtain 

p(l+4p)M'/+(l + 8p)M;+Mi = G, (38) 

where 

G=(l+4p)//; + ^^^//i + -//2- (39) 
P P 

The homogenous part of Eq. (l38l) defines a Gauss hypergeometric function. A change of variable 
\.ok = -\/— 4p, further simplifies the equation 

The two solutions of the homogenous part are 

= K(fc) , (41) 



^2{k) = K Vl-;t2 , (42) 



where K is the complete elliptic function of the first kind defined in Eq. (|2TI) . 

As it stands, the solution Eq. (I42l) requires analytic continuation in the physical range, because 
k is actually imaginary. Using the relation 



fx/T^) = / k( -j^= \ , (43) 



4>2(/:) = JL_ K ( JL_ ) . (44) 



9tK 

we choose to set 

(b^(k\= , ^, , 

At this point, we only need to compute the Wronskian of the solutions in order to be able to 
present a complete solution to our problem. Using the fact that the derivative of the complete 
elliptic function of the first kind is 

dK E{k) K{k) 
dk k{l-k^) k ' 
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Figure 6: Example Feynman diagrams contributing to the two-loop amplitude for heavy flavor 
production in gluon fusion, which need the integral from Fig. \5\without the cut. Thick lines are 
massive as before. 

where E is the complete elliptic function of the second kind defined in Eq. (I22l) . we obtain the 
following Wronskian 

In the above equation we have also used the Legendre relation 

E(yt)K (^Vl -k'^^ +K(yt)E [yi-k^^-K{k)K [^l-k^^ = (47) 
Finally, the solution to Eq. (I38l) is 



Mi(p) = CiK(v^)+C2 



v/4pTT 

1 ri /k(-^)k(v/=4^) K(v/=4p)k(^; 



- / dx 



TtJp \ V4p + 1 V4x+ 1 

xx(4x+l)G(x) . (48) 

The choice of the upper bound above allows to easily find the integration constants C\ and C2. 
After obtaining M2 from Eq. (|34l) . and using the fact that both integrals have to vanish at threshold, 
it turns out that Ci,2 = 0. 

Before we turn to some general conclusions based on the result we just obtained, let us stress 
that in this case the right choice of integration variable has proven to be invaluable. Indeed, had we 
tried to solve the differential equations in our conformal variable x, it would have been impossible 
to identify the differential equations as hypergeometric. 

A quick inspection of Eq. (l48l) shows that it is impossible to express this integral with polylog- 
arithms or a generalization thereof (defined as nested integrals over some elementary functions). 
This statement is in particular true of Goncharov's poly logarithms, and is a simple consequence of 
the definition of the elliptic integrals. In general, the presence of functions with such properties in 
master integrals does not preclude the possibility for their cancellation in the final cross section. 
However, as the results in Section |4] demonstrate, these non-polylogarithmic functions persist also 
in the final result. 

The fact that the total cross section requires a substantially extended set of functions is probably 
not that surprising, because of the never trivial phase space integration. After all, the two masters 
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considered in this section correspond to a three particle cut. However, the analysis we have just 
performed implies that the same functions will be necessary to solve the master integrals for the 
fermionic contributions to the two-loop amplitude for heavy flavor production in gluon fusion. 
Examples of Feynman diagrams that contain the integral from Fig. |5] without the cut are given 
in Fig. [6l In this argument it is crucial that the only place where the presence of a cut plays a 
role are the non-homogenous terms Hi and H2 (and therefore the function G). In consequence, 
the result for the two-loop integral will be given by the very same Eq. (l48l) . This means that it 
will be substantially more difficult to obtain an analytic result in this case than it was in the quark 
annihilation channel discussed in [42]. Finally, let us stress that even though elliptic integrals have 
been encountered previously in Quantum Field Theory (see for instance a recent study in [73]), 
we believe that our observation in this particular problem is an important lesson on complications 
associated with massive QCD processes. Whereas in the massless case ordinary polylogarithms are 
sufficient, with masses present on internal lines one will encounter elliptic functions and possibly 
other yet unexpected objects. 



7 Summary 

In this work we have calculated for the first time the exact analytic expressions for the next-to- 
leading QCD corrections to all inclusive cross sections for heavy quark pair production at hadron 
colliders. The most important phenomenological application is top quark pair production at the 
Tevatron and LHC. 

The same corrections have been computed in approximate form by several groups using more 
traditional methods based on numerical phase space integration. Our calculation utilizes com- 
pletely independent methods where both real and virtual integrations are performed simultaneously 
and in analytic form. We find complete agreement with the earlier calculations; the comparison to 
our exact result shows that the numerical uncertainties in the qq and gq reactions are indeed very 
small. While much larger when compared to the case of the other two reactions, the uncertainty in 
the gg reaction does not exceed the promised accuracy of 1%. 

The motivation for the present calculation was manifold. First, it is a step in the ambitious 
program for the calculation of the NNLO corrections to top quark pair production at the LHC. 
Indeed, UV and coUinear renormalization require the computation of terms subleading in the cor- 
responding regulator, or in the case at hand, beyond d = A. Second, one would like to derive and 
test methods capable of tackling the very demanding NNLO calculation. Third, since the resulting 
complexity is driven by the underlying analytic structures, the current calculation allows one to 
gauge the true level of the complexity of the NNLO calculation. We find for the first time that 
new, a priori unexpected analytic structures do appear even at NLO. What is more important, we 
are able to pinpoint their specific origin. Indeed, at NLO we have to compute a total of 37 master 
integrals and only seven of them have an analytic form differing from what is known from purely 
massless calculations. 



20 



Therefore, the conclusion that we draw from our present work is that the analytic methods 
applied here alow a large degree of automatization and also seem suitable for the calculation of (at 
least) the bulk of the NNLO resuk. 

Finally, the analytic results derived in this paper have important implications to the NLO/NLL 
(and beyond) phenomenology of top quark pair production at the Tevatron and LHC, especially 
for the production close to the partonic threshold. In particular, we demonstrate that although the 
uncertainty in the previous numerical evaluations of the gluon-gluon cross section is within 1% of 
the exact result, the extraction of the so-called "constant" terms that are very important at threshold 
leads to a much larger numerical difference, one that is in fact as large as 7%. 

Our findings for the threshold behavior of the cross section have to be carefully explored, 
especially in light of the fact that the progress and improvements in the top quark pair production 
cross section in the last ten years or so were based on refinements of the behavior of the cross 
section at threshold. A detailed investigation of this effect will be presented elsewhere [34]. 
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A Some well-known results 

For convenience of the reader, we collect in this Appendix known analytic results for different 
contributions to the partonic cross sections, which have not been given explicitly in Section IH 
We begin by presenting the coefficients of the scale dependent logarithms. Notice, that because 
we renormalize the strong coupling constant with n^ + l active flavors, these coefficients are not 
identical to those given in [6]. They can be brought into the same form by decoupling the heavy 
quark, which at this order of perturbation theory corresponds to the simple replacement 

af'+^)(,.2) = a;"'^(^2) (^l + ^^^log(^) +o(a?)j . 

Clearly, decoupling at next-to -leading order only affects the scale dependence. This property does 
not hold at higher orders. In the case of the SU (3) color group, the fij functions read 

/ 484 484 968 4840 1936 \ 

V 243(^+1) ^ 243 (x+ 1)2 ^ 81(;c+l)3 ~ 243(;c+l)4 ^ 243 (x+ 1)5 ) 
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/ 32 32 128 640 256 ■ 

^ V 27(jc+ 1) ^ 27(x+ 1)2 + 27(jc+ 1)3 " 81(x+ 1)4 + 81(;c+ 1)5 J ^ ''''' 

/ 128 128 256 1280 512 

^ V 81(;c+l) ^81(x+l)2^27(x+l)3 ~81(x+l)4^81(;c+l)5 ' ^ '''' 

^ ^ 81(.+ 1)5 



/ 181 383 784 1402 1319 2638 \ 

VT620 ~ 162(x+l) ^ 81(x+l)2 ~ 81(x+l)3 ^ 81(x+l)4 ~ 405(x+l)5 J 
/ 17 _ 61 _ 1 47 28 28 ~ 

^ V27(x+ 1) ~ 54(jc+ 1)2 ~ 27(x+ 1)3 ^ 18(x+ 1)4 " 9{x+ 1)5 ^ 27(x+ 1)6 J ^ ''''' 
/ 2 2 8 4 

^ V 27(jc+ 1) ^ 9(jc+ 1)2 " 27(x+ 1)3 ^ 27(x+ 1)4 J " 
/ 8 8 32 16 \ 

^ V 9(^+^ ^ 3(x+ 1)2 " 9(;c+ 1)3 + 9(x+ 1)4 ) ^^-^'^'"^^ 



2 



181 94 992 693 14743 14743 



2 



360 9(;c+l) 9(;c+l)2 2(x+l)3 36(;c+l)4 90(x+l)5 
/ 173 103 _ 382 1081 _ 269 176 

^ V 18(x+ 1) ~ 6(x+ 1)2 " 9{x+ 1)3 ^ 9(x+ 1)4 ~ 3(JC+ 1)5 ^ 9(x+ 1)6 J ^ ''''' 
/ 7 10 110 155 62 \ 

+ (^^+(^qno2~(^+T)3 + (^TTr~0^^ 

/ 4 23 3 1 ^_ 

^ V (^TTF " 3(jc+ 1)3 + 0^+1)5 "3(jc+ 1)6 J" 

/ 8_ 24 36 28 12 4 ^ 

+ [ ^+ ! + 1)2 (^+1)3 + (^+1)4 (^+i)5 + (^+i)6j^^ 

f ^ 18 32 16 , 

/ 4 12 28 4 12 4 \ 

As mentioned in Section HI neither fgq nor fgq depend on n^. 

To complete our exposition, we also give the Born contributions, keeping the full color depen- 
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dence in concordance with Eqs. (I26l) and (1271) . and setting Tp to its customary value of 1/2 



71 2x(l -x) (jc2 + 4jc+ l) 
A^^^ 3(x+l)5 



71 



'2(x-l)x(x2 + 7x+l) 8x3h(0,jc) 



A^2-l \ ^ \ 3(x+l)5 (x+l)6 

' 2{x - l)x (jc2 + 6x + 1) 2jc (x^ + 8^3 + 6jc2 + Sjc + l) H(0, jc) 



(1A16) 
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